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Abstract 

In these lecture notes some applications of Monte Carlo integration methods in 
Quantum Field Theory - in particular in Quantum Chromodynamics - are introduced 
and discussed. 
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1 Introduction 



The mathematical description of the Standard Model - the theory of elementary par- 
ticle interactions - is based on relativistic Quantum Field Theory (QFT). Relativistic 
QFT is the quantum mechanics of fields defined on the four-dimensional space-time 
continuum. As such it has an infinite number of degrees of freedom - the values of 
field variables in every space-time point. In order to define it, one has to start with the 
quantum theory of a finite number of degrees of freedom: the values of field variables 
in a finite set of discrete points within a finite volume. In most cases the points are 
lattice sites of a regular, hypercubical lattice over a four-dimensional torus. In order 
to define the theory one has to perform the continuum limit and infinite volume limit 
when the spacing of the lattice points goes to zero and the extensions of the torus grow 
to infinity. 

An important simplification from the mathematical point of view is to consider, 
instead of the real time variable, the time to be pure imaginary. In this Euclidean 
space-time the symmetry with respect to Lorentz-transformations becomes equivalent 
to the compact symmetry of four-dimensional rotations and, perhaps even more impor- 
tantly, the quantum mechanical Schrodinger equation is transformed into an equation 
equivalent to the equation describing heat conduction (or e.g. the Brownian motion). 
The consequence is that QFT with imaginary time is equivalent to the (classical) sta- 
tistical physics of the fields. In the Feynman path integral formulation of quantum 
mechanics the exponent in the Boltzmann-factor is the Euclidean lattice action. (Note 
that the "path" in case of the fields is better named as the "history" of the fields in 
the space-time points.) 

The definition of QFT on a Euclidean space-time lattice provides a non-perturbative 
regularization without the infinities which have to be dealt with in perturbation the- 
ory by the renormalization procedure. One can also define perturbation theory on 
the lattice and in this way the lattice gives an alternative regularization for perturba- 
tion theory: the momentum cutoff is implemented by the absence of arbitrarily high 
momentum modes on the lattice. 

The number of discrete points to be considered tends to infinity both in the contin- 
uum limit and infinite volume limit. In order to differentiate between these two infinite 
limits one has to consider the ratio of the effective size of physical excitations to the 
lattice spacing. Obviously, this ratio has to diverge in the continuum limit. In the 
infinite volume limit, on the other hand, the ratio of the size of physical excitations 
to the volume extensions is relevant. In any case, one has to know about the size of 
the physical excitations which is determined by the (bare) parameters in the lattice 
action. In the language of statistical physics, in the continuum limit one has to tune 
the parameters of the lattice action to some fixed point with infinite correlation lengths. 
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If such a fixed point exists, our knowledge in statistical physics suggests universality, 
which means that one can reach the same fixed point (i.e. the same continuum limit) 
with many different lattice actions. 

The most prominent example of relativistic QFT is Quantum Chromodynamics 
(QCD) which is the theory of strong interactions among the six known "flavors" of 
quarks: u-, d-, s-, c- b- and t-quark. QCD is a mathematically closed theory which 
has an unprecedented predictivity: it has only six independent parameters, the quark 
masses. More precisely the parameters of QCD are: m u j ' Aqcd-, n^d/^-QCD, "W A-QCD, 
m c /KQCD, rn b /KQ C D and m t /AQ C D where the A-parameter of QCD Aq C d is an ar- 
bitrary scale parameter of dimension mass. In many applications of QCD only the 
three "light" quarks, the u-, d- and s-quarks are relevant, therefore there are only 
three (small) parameters: m u ^, s / ' ^-qc D ■ All the properties of strong interactions as 
masses, decay widths, scattering cross-sections etc. are, in principle, determined by 
these parameters. 

The somewhat unfortunate circumstance is that, even if in principle determined 
by a very small number of free parameters, it is difficult to tell what are precisely 
the predictions of QCD. The reason is that strong interactions are obviously (at least 
sometimes) strong and therefore calculational methods based on symmetries and on 
perturbation theory only have a limited range of applicability. The only known method 
to evaluate the non-perturbative predictions of QCD theory is lattice QCD. One can 
formulate this in a different way by saying that the validation of QCD as a true theory 
of strong interactions is the task of lattice QCD theorists. 

In this series of (five) lectures on Monte Carlo methods first the different lattice 
formulations of QCD are reviewed (Section [2]). The basic Monte Carlo integration 
methods are introduced in Section [3] and discussed in some detail, including the im- 
portant methods applicable for quark dynamics ( "un-quenching" ) . Section [5] contains 
a selection of some recent developments in order to illustrate recent trends in lattice 
QCD. Finally, the last Section [5] gives a short outlook. 

2 Lattice actions 

The QFT's on the lattice are defined by their Euclidean lattice action. The lattice is 
in most cases a regular, hypercubical one with periodic boundary conditions (torus). 
Lattice elements are the sites (points) and the links connecting neighboring sites. A 
simple case is illustrated by the two-dimensional 4x4 lattice in Figure HJ The lattice 
spacing is usually denoted by a. For the definition of lattice gauge theories like QCD 
the plaquettes consisting of a closed path of four links are important (see Figure [2]). 

The elementary excitations in QCD are the gluons and quarks. The gluons are 
described by a gauge field with elements in the SU(3) color group U Xfl £ SU(3) co i or 
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4x4 periodic lattice 



Figure 1: A two-dimensional periodic 4x4 lattice. 

associated with the links (x — ► x + /t) where £i denotes the unit vector in the direction 
fi (= 1,2,3,4). These are parallel transporters of the color quantum number. The 
corresponding SU(3) Lie algebra element A x ^ can be defined by the relation U Xfl = 
exp(—aA XfJj ) with the lattice spacing a, in order to display the mass dimension of A Xfl . 
The components of A x/1 are introduced by A x/1 = —igA b tl {x)^\f ) , with the Gell-Mann 
matrices A&, (b = 1, ... ,8) and g denoting the bare gauge coupling. The quark fields 
^ and ^ are associated with the lattice sites, as shown in Figure [2j (For notation 
conventions see, in general, the book pp.) 

plaquette in LQCD 

n < n 



Figure 2: The plaquette. 
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2.1 Lattice actions for gluons and quarks 

2.1.1 The plaquette lattice action of the gauge field 

As stated in the introduction, the lattice action for a given theory is not unique. 
There are large varieties of lattice actions in the same universality class realizing in the 
continuum limit the same QFT. For the lattice action of the SU(3) color gauge field 
in QCD the simplest choice is the Wilson plaquette action introduced by Ken Wilson 
in his seminal paper on confinement and lattice QCD [2]. It is based on the definition 
of the field strength F^ u (x) associated with the plaquette variable 

Ux;nv = Ul^U ] x+i> ^U x+ ^U x ^ = exp[-a 2 G^(x)] , (1) 

where 

G flu (x)=F flu (x) + 0(a) (2) 

and 

Fjiv (x) = Af l A u (x) - AfAp (a?) + [A„ (x),A v (x)] (3) 

with the lattice forward derivative defined as A* ip(x) = <p(x + jl) — <p(x). 

As one can easily show, in general, for an SU(N C ) color gauge field we have 

ReTrU x .^ = N c + ^TrF^(x) 2 + 0(a 5 ) (4) 

and therefore the Wilson (plaquette) gauge field action for the SU(N C ) gauge field can 
be defined as 

Sgauge = S g = ^ E # { 1 ~ W Re ^ ( U ^) f 

= -ijE a ^v(^w+^ 5 )' (5) 

c XpuV 

Here we introduced the usual lattice variable for the bare gauge coupling as 

P = ■ ( 6 ) 

An important property of the Wilson action in ([5]) is gauge invariance. This is due 
to the fact that the trace of the product of link variables along any closed path is gauge 
invariant because the gauge transformation of the gauge link variables is 

U' xli = A-\x + jl) U XiX A(x) [A(x) e SU(N C )] . (7) 

The expectation value of some function of link variables 0[U] is given in terms of 
the invariant group (Haar-) measure dU Xfl as 

(O) = i ! Y[dU xlx e W {-S gauge [U}} Op] = f [dU] e - s »~«"M Q[U} , (8) 

J x\x •* 



where the partition function for the gauge field is defined as 

Z = ! Y[dU x ^ e W {-S„[U}} = j [dU] e -*WM . (9) 

This shows that, indeed, in the Euclidean path integral formulation lattice gauge theory 
is equivalent to the statistical physics of gauge fields. 

2.1.2 The Wilson lattice action of fermion fields 

The Dirac equation for fermions can also be similarly discretized as the equations of 
motion for the gauge field. A simple choice is the Wilson action for fermions: 

gWilson = £ J - \Y.^ + ^U X ^ X - \Y$*+* U ** ~ ^* f • ( 10 ) 

x \ fl fi J 



Here ip x , ip x are anticommuting Grassmann variables which have, in general, a Dirac- 
spinor, a color and a flavor index. For a single species ("flavor") of fermions, of course, 
there is just a spinor and a color index. The lattice spacing is set now to unity: a = 1, 
which is often done in the literature, /xo is the bare quark mass in lattice units and 
the Wilson parameter is r ^ 0. The summation in ()10p runs over both positive and 
negative directions: = X)u=±i an d' by definition, we have 7_ M = —7^- The role 
of the Wilson term proportional to r will be discussed below. In (|l(jp the interaction 
of the fermion with a gauge field is introduced by the gauge link variables U xu . Free 
fermions with no interaction correspond to U xa = 1. 

Often used notations are based on redefining the field normalizations according to 

(/i + 4rf/ 2 i> x ip x , Gu Q + 4rf/ 2 ^ => ^ (11) 

and introducing the hopping parameter by 

K =(2 M +8r)- 1 , ^ = _ 8 r) . (12) 

In this way the Wilson action (|l(jp can be rewritten as 



Wilson 



\ ^x^x) - t^fe +A f/ v [r + la]lpx) \ = ^{ipyQyx^x) ■ (13) 
x I M ) xy 

In the second form the Wilson fermion matrix is (without explicit color- and Dirac- 
indices): 

Qyx = 5yx - ft 3y,x+ji U xa (r + 7^) . (14) 

The particle excitations of Wilson lattice fermions can be identified by considering 
the Wilson fermion propagator, which is defined by the inverse of the (free) fermion 
matrix in (|14f) : 

} ^ ^zyQyx — $zx j ^yx — ^y—x — 7^ ^ ~] 6 ^ ^ ^k • (!"-*) 



n 

k 



Here ft = N1N2N3N4 is the number of lattice points and the allowed values of the 
momenta for periodic and antiperiodic boundary conditions, respectively, are 

ap» = k lx = ^-v„, ^ = f^(^ + ^) (^G{0,l,2,...,iV M -l}) . (16) 

Using the notations 

k.. = 2 sin 

2 



k 

kfj, = 2 sin , kfj, = sin k^ , (17) 



the solution of Eq. (fT3j) is given by 

^ = 1 _ rK ( 8 _ p) _ 2 j K1 • k = ! Mo + (r/2)k 2 -i-y-k 

k [l-r K (8-fc 2 )] 2 +4K 2 P [/i + (r/2)*; 2 ] 2 + fc 2 ' 

Particle excitations belong to the poles of the propagator. Considering the Wilson 
fermion propagator in (I18p . it becomes clear why the non-zero value of the Wilson 
parameter r is required, namely, for avoiding additional particle poles at kn = tt 
besides the physical ones at k^ = 0. For r = 0, which corresponds to the naive 
discretization of the Dirac equation, these additional particles emerge and - instead of 
a single fermion flavor - sixteen flavors are described. The 15 extra unphysical particles 
are the consequence of the first order character of the Dirac equation. Introducing a 
non-zero r removes the unphysical fermions from the spectrum in the continuum limit 
(a — » 0) because their masses tend to infinity as a -1 . The price to pay for repairing 
the particle content is, however, rather high because for r / the chiral symmetry is 
broken also for zero fermion mass! 

2.1.3 The Kogut-Susskind staggered lattice action of fermion fields 

As discussed in the previous subsection, the "naive" fermion action without the Wilson 
term (i.e. r = 0) describes 16 fermion "flavors". The naive fermion action is: 

S naive = £ I + 1 £ _ ^ x+ ^ x ] 1 . (19) 

x [ n=l J 

One can perform on this a spin diagonalization by a transformation 

= A x ^ x , V x = (20) 

in such a way that 

Al^A x = a Xfl 1 4 = 1 4 , Gu = 1, 2, 3, 4) . (21) 

One out of four identical components gives the "staggered" fermion action: 



gstaggered = £ J + 1 £ 

x I fl=l 



+A " 1>*+M } ■ ( 22 ) 



The staggered fermion action describes four degenerate flavors with components 
scattered on the points of 2 4 hypercubes. (Note that there are no Dirac spinor indices 
for staggered lattice fermions - only color indices!) Rather remarkably, at zero fermion 
mass /xo = there is a remainder of exact chiral symmetry, namely, U even (l) ®U dd(^) ■ 

2.2 Improved fermion actions 

The freedom of choosing the lattice action in the universality class of the same limiting 
theory in the continuum can be used for: 

• accelerating the convergence to the continuum limit, 

• achieving enhanced symmetries already at non-zero lattice spacings. 

In QCD particularly interesting is the improvement of chiral symmetry at non-zero 
lattice spacings which implies, for instance, simpler renormalization patterns for com- 
posite (e.g. current-) operators. 

The basic tools for constructing improved actions are lattice perturbation theory, 
renormalization group transformations [3] and the local effective theories at non-zero 
cut-off HE]. 

Great effort has been invested recently in constructing improved actions for stag- 
gered quarks (see, for instance, the papers of the MILC Collaboration [6]). In the so 
called Asqtad action the gauge action includes a combination of the plaquette, the 1x2 
rectangle and a bent parallelogram 6-link term. The quark action includes paths up to 
seven links of the form ip y Uy^ x ip x where U yt ^ x is the product of links along the path 
x — ► y. The relative weight of the contributions is such that the flavor symmetry break- 
ing is suppressed and the small momentum behavior is improved. Since one staggered 
quark field describes four "flavors" of fermions (called here "tastes"), for describing a 
single quark flavor in the path integral the fourth root of the fermion matrix is taken 
( "rooting"): 

J [dU d^dip] e~ s ^ = j [dU] e~ s « detQ => J [dU] e~ s * (det Q) 1/4 . (23) 
It is assumed (but debated) that this gives the correct continuum limit. 

2.2.1 Twisted-mass lattice QCD 

A particularly simple way of improving the Wilson-fermion action is the chiral rota- 
tion of the Wilson term in S^ ilson Eq. (QSJ) E]. For two equal mass quark flavors 
(Nf = 2) the unbroken SU(2) subgroup of the SU(2) £g> SU (2) chiral symmetry can be 
partly rotated to axialvector directions. In addition, "automatic" O(a) improvement 
is possible [9]. 



The twisted mass lattice fermion action is: 



tm 



+ a / u cr ^e-^ 5T " 3 ^-^[V', +A ^ At -V'Je-^ 75r3 V'4 • (24) 



cr 



Here u is the twist angle, a^ q the bare quark mass in lattice units and o// cr = (jt, 
4r) < the critical bare quark mass where ijP hysical = o. 

The "twist" can be moved to the mass term by a chiral transformation 

^ i 

Xx = exp(--w7 5 r 3 )^ :2: , x x = i>x exp(- -W75T3) , (25) 

hence the name "twisted mass" . Introducing the quark mass variables 

fi K = afi cr + afx q cos u> = — = amo + 4r , a/i = a\x q sinw , (26) 

the action in ([25!) becomes 

±4 



(xJ/^« + *7573 a/j]Xx) ~ \ E (Xx+a^m^ + 7A*]Xa 



2 aj=±1 



In numerical simulations one starts with this form because it does not contain the 
critical quark mass a\i cr which is a priori unknown and has to be first numerically 
determined. Near maximal twist corresponding to u> = ir/2 it is also convenient to 
introduce till another fermion field by the transformations: 

* x ~ 72 (1 + ' 75T3) Xx ' = Xx^(! + »7573) ■ (28) 

The quark matrix on the x-basis defined in ([27]) is 

1 ±4 

Qxy = S xy (Mk + H57-3 a)u) - - ^, y +/i^[r- + 7 M ] (29) 

/i=±i 

or in a short notation, without the site indices, 

= fi K + Z75T3 o/i + iV + R , (30) 

with 

1 ±4 r ±4 

-^rj/ = ~~ 2 ^ ] fixjy+p.Uyn'Jfj, , ii^j; = — — ^x,y+fiU y ^ . (31) 
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On the ty-basis defined in (I28h we have the quark matrix 

Q W = \( l ~ H5T 3 ) Q ix) (1 - H5T 3 ) =afi + N- i 7 5T 3 (m« + fl) • (32) 

The quark determinant in the path integral over the gauge field is, for instance, 
using the quark mass variables in (I24h : 

det (D cr + a//q cos uj)\D ct + afi q cos a;) + {afi q ) 2 sin 2 uj (33) 

where the single-flavor critical fermion matrix is 

An important feature of the twisted mass formulation is that the fermion matrix 

D cr + a^ g (cos uj + ^75X3 sin uj) (35) 

cannot have zero eigenvalues for non-zero quark mass if uj 7^ 0, tt. There are no spurious 
zero modes and hence no exceptional gauge configurations with anomalously small 
eigenvalues of the fermion matrix. This makes the Monte Carlo simulations at small 
quark- (and pion-) mass easier. 

The consequence of the chiral rotation corresponding to the twist is that the direc- 
tions of vector- and axialvector-symmetries in the SU(2) (g) SU (2) chiral group are also 
rotated. One can achieve conserved axialvector currents but then some of the vector- 
(flavor-) symmetries will be broken. (The twist also induces a breaking of parity.) 
The status and consequences of the chiral symmetry can be deduced from the chiral 
SU(2) ®SU(2) Ward-Takahashi-identities. 

Exactly conserved axialvector currents can be achieved at uj = ^n. In this special 
case the conserved currents are: two axialvector currents (j = 1,2 ) 

A jxl = \{(^x+i,lnl^U x ^ x S) j + (^TmTs^^A+a) 

+ r \^ X+[ ^U X ^ X J - r (^pJj-Ul^ x+ -^j j (36) 

with t\ = T2 and T2 = — T\, and one vector current: 

1 r/- 



*C = ^{{^x + ^^U x ^ x ) + ^ xl ^Ul^ x , 

- ~2 $x+(n5 u *i* 1 f>x) + ~2 {^x^Ul^+fiJ > . (37) 



The invariance of the path integral with respect to the change of variables 

1 1 

V4 = (1 + -^CtVrxTr + -aArx75T r )tpx , 

j 1 1 

i>x = i'xi 1 ~ 2 aVrxTr + 2 a ArxlbTr) (38) 
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implies for an arbitrary function O of field variables the following WT-identities: 

with the backward lattice derivative defined as Aj a <p(x) = tp(x) — cp(x — fi). 

Besides the conserved axialvector currents the important feature of twisted-mass 
Wilson fermions is automatic 0{a) improvement. (0(a) improvement means that in 
the continuum limit a — > the leading deviation from the limiting value behaves 
asymptotically as 0{a 2 ).) As it has been shown by Prezzotti and Rossi [9j, for the 
(untwisted) Wilson fermion action we have 

(°>E) = \ [(0)(r,m q) + (0) { - r , mq) ] OC (0)%% + 0(a 2 ) . (40) 

This is averaging over opposite sign Wilson parameters: "Wilson average". 

In twisted mass lattice QCD (tmLQCD) changing the sign of r is equivalent to 
shifting the twist angle by ir. In the special case of u> = ^tt this is equivalent to u> — > —ui, 
therefore expectation values even in oj are "automatically" O{o) improved, without any 
averaging. Automatically 0{a) improved physical quantities are, for instance: 

• the energy eigenvalues, hence the masses; 

• on-shell matrix elements at zero spatial momenta; 

• matrix elements of operators with parity equal to the product of the parities of 
the external states. 



2.2.2 Domain wall lattice fermions 

The chiral symmetry of massless fermions can be realized at non-zero lattice spacing 
by introducing a fifth "extra dimension" [TT| \\2\ . In the fifth direction there is a 
"defect": either the mass term changes sign [10J or there are "walls" at the two ends 
[12j . In this case there are chiral fermion solutions which are exponentially localized 
in the fifth dimension near these defects. The gauge field remains four-dimensional 
(independent on the fifth dimension) . In the limit of infinitely large fifth dimension the 
positive and negative chirality solutions (at opposite walls or at opposite sign changes 
on a torus) have zero overlap with each other and the chiral symmetry becomes exact. 
The domain wall fermion action can be written (with 1 < s < N s ) as 

S F = ^f«(D f ) ISilV f A - (41) 

s,s' 
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where in an s-block form 

( a + D -aP L 

-oP R a + D -aP L 

-aP R a + D -aP L 

D F 









m f P R \ 





(42) 



... -oP R a + D -aP L 

\ m f P L ... -aP R a + D J 

The chiral projectors are denoted, as usual, by Pr % l = ^(1 ± Is), the quark mass in 
lattice units is rrif, the ratio of lattice spacings is a = a/a s and the four-dimensional 
Wilson-Dirac matrix with negative mass (0 > —niQ > —2) is, for r = 1, 

4 



D xx > 



(4 - m Q )8 z 



(43) 



The hermitian fermion matrix corresponding to Dp in (J42J) is useful, for instance, 
in Monte Carlo simulations. It can be constructed as follows: since with an s-reflection 

(Ik)as> = $N s +i-s,si we have 



Dp = R blb D^ F R hlb 
the hermitian fermion matrix can be defined as 

Dp = R 5T oD F = D F 



(44) 



(45) 



The chiral symmetry is broken by a non-zero overlap of the opposite chirality wave 
functions, which tends to zero in the limit of an infinite extension of the fifth dimension: 
N s — > oo. Enhanced symmetry breaking occurs if the four-dimensional Wilson fermion 
matrix D has small eigenvalues. 

2.2.3 Neuberger overlap fermions 

Another possibility to achieve chiral symmetry of the lattice fermion action, which in 
fact can be related to domain wall lattice fermions, is the Neuberger ( overlap- ) fermion 
action. 

Let us rewrite the (free) Wilson fermion action for r = 1 and /j,q = amo as 



S. 



Wilson 



D W = £ 



^(V M + V*)-^V*V M 



(46) 
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where the lattice derivatives are now denoted by 

V„ = a- 1 A* , V; = a" 1 A* . (47) 

The Neuberger lattice fermion operator with zero mass is defined as 

D N = - [1- A — ) , A=l-aD w . (48) 



The inverse square-root here can be realized by polynomial or rational approximations. 
Note that A is proportional to the Wilson fermion matrix with bare mass — a . 

An important property of the Neuberger operator is that V = 1 — aD^ is 
unitary: V'V = 1. As a consequence, the spectrum of Dn = a~ l (l — V) is on a circle 
going through the origin. In addition, the Neuberger operator satisfies the Ginsparg- 
Wilson relation 

-y 5 D N + -Dtv75 = aD N -f 5 D N . (49) 
This is equivalent to the condition as introduced by Ginsparg and Wilson (GW) [13] 

l5 D- 1 + ITS = 2aR l5 . (50) 

The GW-relation is the optimal approximation to chiral symmetry which can be real- 
ized by a lattice fermion operator for a — > 0. R in (I50p is, in general, a local operator. 
For the Neuberger operator D = Dn we have R = \. 

The lattice chiral symmetry satisfied by a GW-lattice fermion can be explicitely 
displayed by appropriately defined chiral transformations [2] . It can be shown that 

^ = 75(l-^)^, # = ?(l-V) 7 5 (51) 

is an exact chiral symmetry for any lattice spacing a if the GW-relation is satisfied. 
Lattice actions satisfying the GW-relation are: 

• the fixed point action, which is the fixed point of some renormalization group 
transformation |15j : 

• the Neuberger action Dn in 



• the effective (four-dimensional) action of the light fermion field of the domain 
wall fermion [16J. 

Note: the inverse of the effective Dirac operator of the light fermion field of the 
domain wall fermion is equivalent to the inverse of the truncated overlap Dirac operator 
(except for a local contact term) . Using GW-fermions one can prove the index theorem 
about topological charge [17] and introduce the ^-parameter in QCD, etc. 

Having lattice actions with exact chiral symmetry at non-zero lattice spacing is a 
great achievement. Although it is expected that (spontaneously broken) chiral sym- 
metry is restored in the continuum limit also for simple lattice formulations with, for 
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instance, Wilson fermions, the explicit breaking of chiral symmetry for non-zero lattice 
spacings makes the renormalization of composite operators more involved and in prac- 
tice also much more cumbersome because of the extended mixing pattern. The chiral 
symmetry restricts the mixing to be simpler and more tractable. 

The difficulty of defining chiral symmetric lattice actions is emphasized by the 
Nielsen-Ninomiya theorem [TS]. This theorem states that there is no (free) lattice 
fermion action which can be written in the form 



• its Fourier-transform is Dip) = i^^Pfj, + 0(ap 2 ) for p <C 7r/a, 

• D(p) is invertible for p ^ (i.e. there are no massless fermion doubler poles), 

• 75D + D75 = (chiral symmetry). 

GW-fermions circumvent the Nielsen-Ninomiya theorem by relaxing the last condi- 
tion: instead of exact anticommutativity only a weaker condition, namely the Ginsparg- 
Wilson relation in (|49p . is satisfied. Correspondingly, the chiral transformation is mod- 
ified: the simple continuum transformation is generalized to (|5ip . 

The important question is whether the locality of the action is ensured for GW- 
fermions. In case of the Neuberger (overlap) action locality can be proven if the gauge 
field is smooth enough, namely if every plaquette value is close to unity [19]. Because 
of the importance of locality such gauge fields are sometimes called "admissible". Of 
course, usual lattice actions typically admit any plaquette value and therefore in the 
path integral "inadmissible" configurations also occur. In fact, in actual simulations 
there are always plaquettes with small values. It is an open question whether this turns 
out to be a problem in the continuum limit. In any case, the lattice spacing has to be 
small enough in order to avoid the "Aoki phase" with lots of small eigenvalues of Dy/. 
The small eigenvalues make non-local and the "residual mass" breaking the chiral 
symmetry of domain wall fermions large |20j. 



The goal of numerical simulations in Quantum Field Theories (QFT's) is to estimate 
the expectation value of some functions A[tp] of the field variables generically denoted 
by [(f \ = {y>xa}- In terms of path integrals this is given as 




(52) 



.''.V 



and which would simultaneously satisfy the following conditions: 
• D(x) is local (bounded for large x by e _7 ' x '), 



3 Monte Carlo integration methods 





(53) 



14 



S[ip] is the lattice action, which is assumed to be a real function of the field variables. 
(To begin with, we only consider bosonic path integrals.) 

A typical lattice action contains a summation over the lattice sites. Since the 
number of lattice points f2 is large, there are many integration variables. However, since 
(I53p corresponds to a statistical system with a large number of degrees of freedom, in 
the path integral only a small vicinity of the minimum of the "free energy" density will 
substantially contribute. A suitable mathematical method to treat with such situations 
is Monte Carlo integration. (For a recent review of Monte Carlo integration in QFT's 
see Ref. [21].) 



3.1 Monte Carlo integration 
3.1.1 Simple Monte Carlo integration 

Let us consider a continuous real function f(X) of a continuous random variable X 
having probability distribution px{s) and hence the expectation value 

(f(X)) = J dsf(s)p x (s) . (54) 

Using px(s) to obtain N outcomes of X (Xx, X2, ■ ■ ■ , Xpj), the random variables Yj = 
f(Xj) give 



1 - f 
lim - Y,Yj = (Y) = (f(X)) = / dsf(s)p x (s) . (55) 



In a short notation: 

N 



5=1 

For large N, the central limit theorem tells us that the error in approximating 



(f(X)) is given by the variance V[f(X)} as y/V[f(X)]/N. The Monte Carlo estimate 
of the variance is: 



V[Y] = ((Y - (Y)) 2 ) « (/ - ff = p-f. (57) 

Generalizing this to several {D) integration variables one obtains the following formulas 
for simple Monte Carlo integration: 

^d D xp(x)f(x)^J± (^j^J ■ (58) 

Here, according to the notation introduced in (f56l) . 

I v I v 

/^E/(^)> / 2 -^E/(^) 2 - (59) 

1=1 1=1 
The points a?i, £2, ■ ■ ■ , %n have to be chosen independently and randomly with proba- 
bility distribution p(x) in the D-dimensional volume V. 
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3.1.2 Importance sampling 

Simple Monte Carlo integration works best for flat functions but is problematic if the 
integrand is sharply peaked or rapidly oscillating. Therefore, a good procedure is to 
apply importance sampling: find a positive function g(x) with integral norm unity 
(J dxg(x) = 1) such that h(x) = f(x)/g(x) is as close as possible to a constant and 
then calculate 

f 'dxf(x) = f dxg(x)h(x) «^7^£M^) , (60) 

J a J a 

where the points Xj are chosen with probability density g{x) and we used simple Monte 
Carlo integration with a constant probability in an interval: 

• (61) 

The prerequisite is, of course, that one can find an appropriate g(x) such that on can 
generate points with it. 

How can one generate the desired (in general, multi-dimensional) probability dis- 
tributions? One possibility for lower-dimensional integrals is the rejection method. 
This is based on the observation that sampling with px(x), for instance, in an interval 
x € [b, a] is equivalent to choose a random point uniformly in two dimensions and reject 
it unless it is in the area under the curve px(x). For high-dimensional distributions 
this becomes cumbersome. Multi-dimensional integrals can be handled by exploiting 
Markov processes. 



/ 

J a 



3.1.3 Markov chains 

A Markov process (or "Markov chain") is a sequence of states which are generated with 
transition probabilities from a given state to the next one. The transition probability 
is assumed to depend only on the current state of the system and not on any previous 
state. For simplicity, for discrete states s%, S2, • • • , sr the transition probability can 
be denoted by pij. The matrix P with elements pij is called transition matrix (or 
Markov-matrix) . 

The mathematical properties of Markov chains are extensively covered in the liter- 
ature. For a comprehensive collection of features relevant in Monte Carlo integration 
of QFT's see Ref. |21j . Let us mention here just a few of them: 

• The product of two Markov matrices P1P2 is again a Markov matrix. 

• Every eigenvalue of a Markov matrix satisfies |A| < 1. 

• Every Markov matrix has at least one eigenvalue A = 1. 
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A very important statement is given by the fundamental limit theorem for (irre- 
ducible, aperiodic) Markov chains: they have a unique stationary distribution satisfy- 
ing w T = w T P which is identical to the limiting distribution Wj = lim r j_» 00 p^' . 

An important concept is the autocorrelation in Markov chains. Since the state of 
the system depends on the previous state, the consecutive states are not uncorrelated. 
To reach a more or less uncorrelated distribution from some initial one, in general, 
several steps have to be performed. The degree of correlation among the subsequent 
states can be characterized by the autocorrelation function which is defined for some 
observable Oi as 

p{t) = ((O t O t+t ) - (O t ) 2 ) /((Of) - (O t ) 2 ) . (62) 

Obviously, decreasing autocorrelations decrease the Monte Carlo error for a given 
length of the Markov chain. 



3.2 Updating 

The aim in Monte Carlo simulations of QFT's is to calculate the expectation values of 
some functions of field variables as given in (|53p . The Monte Carlo integration is based 
on importance sampling. The required distribution of field configurations according to 
the Boltzmann factor e~ s ^ ( "canonical distribution") is generated by a Markov chain 
by exploiting the fundamental limit theorem discussed in Section 13.1.31 

Let us denote the configuration sequence generated in the Markov chain by {[</?«], 1 < 
n < N}. In this field configuration sample the expectation values are approximated by 
the sample average: 

A =jfY, A M (A)- (63) 

n=l 

The Markov process of generating one field configuration after the other is generally 
called updating. Let us denote the transition probability from a configration to the next 
one [ip] —* [ip'] by P ([(p'] <— [<p])- In order to generate the canonical distribution e~ s ^ 
a sufficient condition is 

P ([p'] - [ip]) e~ s M = P ( M <_ [</]) e -m . (64 ) 
This condition is called detailed balance. 



3.2.1 Metropolis algorithm 

The "ancestor" of updating processes for bosonic systems is the Metropolis algorithm 
|22j . For a system with N possible configurations the transition probability for [<p r ] ^ 
[ip] is defined by 

P(W\ - M) = N min \ 1, e — m \ . (65) 
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This transition matrix can be realized by the following numerical procedure: 

i. ) choose first a trial configuration randomly from J\f configurations and 

ii. ) accept it as the next configuration in any case if the Boltzmann factor 
is increased (the action is decreased). If the Boltzmann factor is decreased 
(the action is increased), then accept the change with probability equal to 
the ratio of the Boltzmann factors. 

The accept-reject step can be implemented by comparing the ratio of the Boltzmann 
factors to a pseudo-random number between and 1. One can see by inspection that 
the above transition probability distribution satisfies the detailed balance condition 
(|64p . hence it creates the desired canonical distribution of configurations. 



3.2.2 Fermions in Monte Carlo simulations 

The lattice action for QFT's with fermions, for instance like QCD, has the generic form 

S[U,iP,^]=S g [U} + S q [U,tP,^] , (66) 

where S g is the bosonic part, in QCD the color gauge field part, and S q is describing 
the fermion fields and their interaction with the bosonic fields. S q is assumed to be 
quadratic in the Grassmann- variables of the fermion fields: 

Sq = ^(^yQyx^x) • (67) 
xy 

The expectation values have the general form 

After performing the Grassmann integration one obtains 

tyAMx* ■ ■ ■ i>v3 Xn F[U\) = Z- 1 J [dU]e- s *M det Q[U] F[U] 

■ E e E£"2 QMwQP]*** ■ ■ ■ Win • (69) 

Zi—Zn 

Here Q[U] _1 is an (external) quark propagator and det Q[U] generates the virtual 
quark loops. 

Since taking into account the fermion determinant det Q [U] in the path integral 
over the bosonic (gauge-) fields is a very demanding computational task, in a crud ap- 
proximation one sometimes simply omits it. This is called "quenched approximation": 
det Q[U] 1. Experience in QCD shows that the results in the quenched approxima- 
tion are often qualitatively reasonable, nevertheless the error caused by omitting the 
closed virtual fermion loops is uncontrollable and implies the presence of unphysical 
"ghost" contributions. 
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3.2.3 Dynamical fermions: "unquenching" 

In the early days of lattice QCD simulations quite often the quenched approximation 
was taken. This is, however, on the long run not acceptable, the obtained results do 
not represent a numerical solution of QCD. More recently - due to some impressive 
developments in the available computer power and in our algorithmic skills - the true 
dynamical simulation of quarks became feasible. 

The basic difficulty in "unquenching" is that the fermion determinant is a non-local 
function of the bosonic fields and therefore it is a great challenge for computations. 
For solving this problem a useful tool is the pseudofermion representation |23j : 

det (QtQ) oc j [##+]exp j- ^(^[Q+Q]^) j . (70) 

In case of, for instance, Wilson quarks the quark determinant satisfies 

Q ] = 75Q75 => det Q f = det Q , (71) 

therefore Eq. (|7ip describes the quark determinant of two degenerate quark flavors. 

In the popular Hybrid Monte Carlo (HMC) algorithm [23] the representation (fT0|) is 
implemented in the updating by using molecular dynamics equations (see Section 13.31) . 
For single quark flavors HMC is not applicable. One can, however, use Polynomial 
Hybrid Monte Carlo (PHMC) [251 HI] (see Section E3D or Rational Hybrid Monte 
Carlo (RHMC) [27]. 



3.3 Hybrid Monte Carlo 
3.3.1 HMC for gauge fields 

The basic idea of HMC is to employ molecular dynamics (MD) equations in order to 
collectively move the field configuration in the whole lattice volume. Since discretized 
molecular dynamics equations are used, the lattice action (analogous to the energy in 
molecular dynamics) is not conserved along MD -trajectories, therefore at the end of a 
trajectory a Metropolis accept-reject step has to be implemented. In this subsection 
HMC will be introduced in the important case of lattice gauge fields, specifically SU (3) 
(color) gauge field. 

The equations of motion are derived from a Hamiltonian which is defined for the 
colour gauge field U X)fi G SU(3) as 

H[P,U\ = l -Y.P* M + S M > ( 72 ) 
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where S g [U] is the gauge field action and the real variables P Xf ij, j = 1, . ■ ■ ,8 are called 
conjugate momenta. They are the expansion coefficients of the Lie algebra element 

j 

It is assumed that the conjugate momenta have a Gaussian distribution: 

P XM oc exp | - Pl M \ = P M [P] . (74) 

The expectation value of some function F[U] is defined as 

= f[dP][dU]exp(-H[P,U])F[U] 
{ ' J[dP}[dU]exp{-H[P,U}) ' 1 ' 

By a proper choice of the discretized trajectories one can achieve that the tran- 
sition probability from a configuration to the next satisfies detailed balance (see next 
subsection). Therefore, the correct canonical distribution is reproduced. 

The Hamiltonian equations of motion are: 

^ = -D Xflj S g [U] , ^ = iP^Us* , (76) 

where the derivative with respect to the gauge field is defined, in general, as 

d 



/(e iQA ' E4,„) . (77) 

a=0 



3.3.2 Detailed balance 

In order to prove that HMC reproduces the correct canonical distribution of (gauge) 
fields it is sufficient to prove the detailed balance condition (|64p for the transition 
probabilities realized by the MD-trajectories. 

The discretized trajectories Tjj provide the following transition probability distri- 
bution at the end of the trajectory: 

P H ([P', U'} <- [P, U\) = 5 {[P', U'\ - T H [P, U)) . (78) 

Let us assume that the trajectories satisfy reversibility: 

P H [[P', U'] - [P, U}) = P H [[-P, U] - [-P>, U'}) . (79) 

The Metropolis acceptance step is described by the well known probability distribution 
p A Qp', U>] _ [ P) [/]) = min |l, e -H[P',U']+H[P,U] j (8Q) 
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The total transition probability is then 

P ([i/'] _ [[/]) = | [ d PdP']P A ([P', 17'] - [P, 17]) Ph ([P', U'] - [P, E7]) Pm[P] • 

(81) 

Using the relation 

e -^] m in{l, e -^'^']+^^} =e-^ P '^min{l,e- H ^ +H ^ P '^} , (82) 
one shows 

= e-^'^P^a-P^l-hP',^]) . (83) 
Therefore, due to reversibility, we have for the canonical distribution 

W C [U] ocexp{-S g [U]} (84) 

the relation 

W C [U] J [dP dP']P A ([P', [/'] <- [P, U]) P H ([P', U'} <- [P, 17]) P M [P] (85) 
= W C [U'} j [dPdP']P A ([-P, E/] <- [-P', U']) P H ([-P, U] - [-P', [/']) Pm[-P'] • 
Taking into account that 

[dPdP'\ = [d(-P)d(-P')} , (86) 
this is just the detailed balance condition we wanted to prove. 

3.3.3 Leapfrog trajectories 

The proof of detailed balance for HMC in the previous subsection has been based on the 
assumption that the discretized MD-trajectories are reversible. The classical example 
is a leapfrog trajectory which is defined as follows. 

First we update the conjugate momente with a step size At = \o~t. This is followed 
by (n — 1) update steps with At = ^St both for the gauge variables and for the 
momentum variables, alternating with each other. Finally, the gauge variables are 
updated with At = 5t and the momentum variables with At = \o~t. 

The explicit formulae for these steps are: 

= Pxfji'j — Dx^j S g [U] At 

= expj^iAj-P^ArJl/^. (87) 
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One can easily prove that the reversibility condition (I79p is satisfied. 

The single steps in the leapfrog trajectory cause a discretization error of the order 
<5t 3 . Therefore, the action for the final configuration is expected to differ from the 
initial configuration by an error of order 5t 2 . 

In the second equation of (I87p we need, in each step on a trajectory for each link, 
the evaluation of the exponential of an element of the gauge group algebra A. It is 
desirable to minimize the cost of this, but at the same time the calculation has to be 
precise enough for not loosing reversibility. Since one can show that 

A 3 = A+ Qt>^ 3 ) I , (88) 

any analytic function f(A) can be written as 

f(A) = a 2 A 2 + ai A + a I . (89) 

For the exponential function f(A) = exp(A) the coefficients ao,i.2 can be practically 
calculated by recursion relations based on the Taylor expansion of exp(A). 

3.3.4 HMC for QCD 

Besides the color gauge field dealt with in the previous subsections, in QCD one has 
to introduce the quarks, too. Let us consider here two equal mass quarks, in order to 
be able to replace the fermionic quark fields by bosonic pseudofermion fields according 
to (jTOj) . (Single quark flavors will be considered in the next Section l3~H ) 

Let us note that the pseudofermion field in (170)) is an auxiliary complex scalar field 
4>qxac having the same number of components as the fermion field i/jqxac- (The indices 
in QCD are: q for the quark flavors, x for lattice sites, a for the Dirac spinor index 
and c for color.) According to (I70p the fermion determinant induces an effective action 
for the gauge field which can be written as 

SeffP] = J2(ct>t{Q{U] + Q[U]}^<j>*) • (90) 

xy 

In the MD-trajectories of the previous subsections S e ff[U] has to be added to the pure 
gauge action: 

Sg[U]=>Sg[U]+S e ff[V] • (91) 

3.4 Polynomial Hybrid Monte Carlo 

Here we discuss the PHMC algorithm [25] with multi-step stochastic corrections |26j. 
This update algorithm is applicable for any number of quark flavors, provided that the 
fermion determinant is positive, which is the case for positive quark mass. For negative 
quark masses there is a sign problem, which will not be discussed here. 
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For Nf = 1,2,... degenerate quarks one uses 

at r +i N f/ 2 r ~o i N r/ 2 1 
|det(Q)|"' = {det(Qt Q) \ = {det(Q 2 )} ~ , (92) 

L J L J detP n (Q 2 ) 

where the Hermitian fermion matrix is Q = 75 Q an d the polynomial P n satisfies 

lim P n (s) = x~ N f /2 (93) 

n— >oo 

in an interval [e, A] covering the spectrum of Q^Q = Q 2 . 

The effective gauge action representing the fermions in the path integral is now 

xy 

Sometimes it is more effective to simulate several fractional quark flavors: 



N f /2 

detQ zl 



det Q' 



N f /(2n B )' 



(95) 



which can be called determinant break-up. In this case we need a polynomial approxi- 
mation 

P n (x) ~ x~ a (96) 

with 

a = —L (97) 
2n B 

and positive integer ns- The effective gauge action with determinant break-up has 
then multiple pseudofermion fields: 

n B 

SeffP] = J2 Y.^tyPn(Q 2 )y^ kx ) • (98) 
k=l xy 

Since polynomial approximations with a finite n cannot be exact, one has to correct 
for the committed error. One can show that for small fermion masses in lattice units 
am <C 1 the (typical) smallest eigenvalue of Q 2 behaves as (am) 2 and for a fixed quality 
of approximation within the interval [e, A] the degree of the polynomial is growing as 

n oc \fe oc (am)" . (99) 

This would require in realistic simulations very high degree polynomials with n > 10 3 - 
10 4 . The way out is to perform stochastic corrections during the updating process 



This goes as follows: for improving the approximation a second polynomial is in- 
troduced according to 

P 1 (x)P 2 (x) ~ x~ a , x £ [e, A] . (100) 
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The first polynomial P±(x) gives a crude approximation 



Pi(x)~ x - a . (101) 

The second polynomial P 2 (x) gives a good approximation according to 

P 2 (x) ~ [x a P 1 (x)\- 1 . (102) 

(This can also be extended to a multi-step approximation [26J.) 

During the updating process Pi is realized by PHMC updates [25], whereas P 2 is 
taken into account stochastically by a noisy correction step. This goes as follows: one 
generates a Gaussian random vector with distribution 

p -r?t P2(Q[U] 2 )v 

(103) 



f[drj]e-v f P2(Q[U]^ v 
and accepts the change [U] —* [[/'] with probability 

min{l,A(7 ? ,[[/'] <- [17])} , (104) 

where 

Afa, [[/'] <- [C/]) = exp {-^(Qft/'] 2 )?? + rytp^Q^] 2 )^} . (105) 

It can be shown that this update procedure satisfies the detailed balance condition. 

The Gaussian noise vector rj can be obtained from rf distributed according to the 
simple Gaussian distribution 



e 



-77't rj' 



(106) 



by setting it equal to 

V = P2(Q[U] 2 )-^v' ■ (107) 

In order to obtain the inverse square root on the right hand side one can proceed with 
a polynomial approximation 

P 2 (x) ~ P 2 (x)-^ , x G [e, A] . (108) 

The interval [e, A] can be chosen differently from the approximation interval [e, A] for 
P 2 , usually with e < e. 

The polynomial approximation with P 2 can only become exact in the limit when 
the degree n 2 of P 2 is infinite. Instead of investigating the dependence of expectation 
values on n 2 by performing several simulations and extrapolating to n 2 — > oo, one fixes 
n 2 to some high value and performs another correction in the expectation values by still 
finer polynomials. This is done by reweighting the configurations. This measurement 
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correction is based on a further polynomial approximation P' with degree n' which 
satisfies 

lim P 1 {x)P 2 {x)P'(x) = x~ a , x £ [e',X] . (109) 

n'— >oo 

The interval [e', A] can be chosen such that e' = 0, A = A max , where X max is an absolute 
upper bound of the eigenvalues of Q 2 . 

In practice it is more effective to take e' > and determine the eigenvalues below 
e' and the corresponding correction factors exactly. For the evaluation of P' one can 
use recursive relations, which can be stopped by observing the required precision of 
the result. 

After reweighting the expectation value of a quantity A is given by 

(exp^ni-PW)]^ ' 

where r\ is a simple Gaussian noise. Here (. . .)jjn denotes an expectation value on the 
gauge field sequence, which is obtained in the two-step process described before, and 
on a sequence of independent rj 's of arbitrary length. 

In most practical applications of PHMC with stochastic correction the second step 
(or the last step if multiple correction is applied) of the polynomial approximation can 
be chosen precise enough such that the deviation from the exact results is negligible 
compared to the statistical errors. In such cases the reweighting is not necessary. How- 
ever, for very small fermion masses reweighting may become a more effective possibility 
than to choose very high order polynomials for a good enough approximation. 

A positive aspect of reweighting is related to the change of the topological charge 
of the gauge configurations. Such changes occur through configurations with zero 
eigenvalues of the fermion determinant where the molecular dynamical force becomes 
infinite. This implies an infinite barrier for changing the topological charge which 
may completely suppress transitions between the topological sectors. This problem is 
substantially weakened by PHMC algorithms because the polynomial approximations 
do not reproduce the singularity of the inverse fermion determinant (i.e. the zero 
of the determinant) [28J. In this way the gauge configuration can tunnel between 
topological sectors. The more frequent occurrence of the configurations near the zeros 
of the fermion determinant is corrected by the reweighting. 

3.4.1 PHMC and twisted mass 

Until now we tacitly assumed that we use ordinary ( "untwisted" ) fermions. In case of 
twisted mass lattice QCD the numerical simulation of light quarks is, in fact, easier, 
because the quark determinant of a degenerate quark doublet becomes, according to 
Eq. J33D, 

det(Q 2 + ^) (HI) 
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where fi s = fjL q smoj with \i q the quark mass in lattice units and u> the twist angle. 
The polynomials P\, ni (x) and P2,n 2 ( x ) now satisfy 

lim P Xm (x) P 2 . n2 (x) = (x + nl)- N 'l* , x G [e, A] . (112) 

ri2— >oo 

In case of to ~ \ the polynomial approximations have lower orders and the updating is 
faster because of the absence of exceptional configurations with very small eigenvalues, 
due to the presence of the lower limit [j? s . (Note that the very small eigenvalues are 
often originating from topological defects at the cutoff scale, which are unphysical lattice 
artifacts going away in the continuum limit.) 

4 Some recent developments 

In spite of substantial algorithmic developments, lattice QCD simulations near the 
small (physical) quark masses still need rather high computer power: we need Tflops! 
An example for a demanding Monte Carlo simulation (in the near future) is: O = 
50 3 • 100 = 1.25 • 10 7 and am q = 0.005. This is equivalent, for instance, at a = 0.1 fm 
to m q = lOMeV, L = 5 fin, ~ 200 MeV. 

The smallness of the u-, d- and s-quark masses implies that the numerical simulation 
(with dynamical quarks) is a great challenge for computations. There are a number of 
large international collaborations working on this problem over the world: 

• USA: MILC, RBC, ... Collaboration; 

• Japan: CP-PACS, JLQCD, ... Collaboration; 

• Europe: UKQCD, Alpha, QCDSF, ETM ... Collaboration. 

It would be rather difficult to give a review of all the interesting results achieved 
over the last years. Here I shall only give a very limited and personal collection of some 
of the problems and results. 

4.1 The light pseudoscalar boson sector 
4.1.1 Gasser-Leutwyler coefficients 

The physical consequence of the smallness of three quark masses is the existence of 
eight light pseudo-Goldstone bosons: tt, K, n. In the low-energy pseudo-Goldstone 
boson sector there is an SU(3) ® SU(3) chiral flavour symmetry and the dynamics 
can be described by Chiral Perturbation Theory (ChPT) [291 130]- In an expansion in 
powers of momenta and light quark masses several low energy constants - the Gasser- 
Leutwyler constants - appear which parameterize the strength of interactions in the 
low energy chiral Lagrangian. 
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An eminent task for Monte Carlo simulations in Lattice-QCD is to describe the 
pseudo-Goldstone boson sector. The Gasser-Leutwyler constants are free parameters 
which can be constrained by analyzing experimental data. In the framework of lattice 
regularization they can be determined from first principles by numerical simulations. 
In numerical simulations, besides the possibility of changing momenta, one can also 
change the masses of the quarks. 

ChPT can be extended by changing the valence quark masses in quark propagators 
independently from the sea quark masses in virtual quark loops. In this way one arrives 
at Partially Quenched Chiral Perturbation Theory (PQChPT) [31] (see Section [4. 1.3ft . 

4.1.2 E(uropean) T(wisted) M(ass) Collaboration 

This collaboration consists of about 30 physicists from 7 countries: 

1. Cyprus: University of Cyprus, 

2. Prance: Universite de Paris Orsay, 

3. Germany: DESY, Universitat Minister, TU Miinchen, 

4. Italy: Universita di Roma 1,11,111, INFN, ECT*, 

5. Spain: Universidad Valencia, 

6. Switzerland: ETH Zurich, 

7. United Kingdom: University of Liverpool. 

In a recent paper (first of a series) numerical Monte Carlo simulations on "Dynamical 
Twisted Mass Fermions with Light Quarks" are reported [32J. 

As examples of the results, Chiral Perturbation Theory (ChPT) fits of the pseudo- 
scalar- (pion-) mass (in Figure[3|) and pseudoscalar- (pion-) decay constant (in Figured]) 
are shown. It is remarkable that the precision on Z3 4 is much higher than obtained 
by any previous experimental determination. However: this is with only Nj = 2 
degenerate dynamical quarks (u- and ci-quark) and no continuum limit extrapolation 
is yet performed (it is comming soon). 

4.1.3 Ratio tests of PQChPT 

Taking ratios at fixed gauge coupling (/?) is advantageous because the Z-factors of 
mutiplicative renormalization cancel (for instance, in m q and f n ). Also: some types of 
lattice artifacts may cancel. 

In case of simulations with Wilson-type lattice actions, by taking into account 
lattice artifacts in the Chiral Lagrangian, one can reach the continuum limit faster. 
This approach is based on the effective continuum theory introduced by Symanzik [3]: 
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Figure 3: Left: (am^) 2 as a function of the twisted mass a^i; right: 
(am T ) 2 /(a/i) versus a\i (by the ETM Collaboration). The finite volume 
ChPT-fit is shown, together with the infinite volume limit (dashed line): 
l 3 = 3.65(12). 

cutoff effects (in the lattice regularized theory) can be described by C(o,a 2 ,.--) terms 
in a local effective Lagrangian. 

This idea can be applied to low energy LQCD [331 E3] - In case of the Wilson quark 
action the leading 0{a) effects have a simple chiral transformation property, identical 
to those of the quark masses. At leading order of ChPT, besides the quark mass 
variable %, an additional 0(a) parameter p appears: 

xs —- ps ir v x)- im 

At next to leading order (NLO): the Gasser-Leutwyler constants L±, . . . , Ls are doubled 
by the (bare parameter dependent) coefficients Wx,...,W% describing 0(a) effects. 
(Extension to 0(a 2 ) is possible.) 

Variables to be used in ratio tests of PQChPT (the index V always stands for 
"valence" quarks which are "quenched", S for dynamical "sea" quarks): 

(i) 

i -_m qV xv _ PS _ m q s Xs M .,s 

£ = = , VS = , 0- i= — -jrr = . 114) 

m qS XS XS m ( 7 XR 

qo 

For the pion decay constants the appropriate ratios are: 

Rfvv^^, Rfvs^^, RRf^jtij-, (115) 

JSS JSS JVVJSS 
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Figure 4: ChPT fits to af n versus afi (by the ETM Collaboration). Left: 
the point with largest afi left out (the dashed line is the infinite volume 
limit); right: compared to finite volume fit to every point. The fit gives: 
a = 0.087(1) fm, a" 1 = 2264(26) MeV), l 4 = 4.52(06). 



and for the pion mass-squares (dividing by the leading order behaviour): 
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are appropriate. 

Examples of the NLO formulas are [341 135j : for N s degenerate sea quarks 
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and similarly for Rn .... 

In the above formulas Lgk denote Gasser-Leutwyler constants renormalized at the 
scale foy/XS- They are related to L k defined at the scale /o and L' k defined at the 
generic scale \i according to 



Lsk = L k -c k log(x< 



L' k -c k logM X s) 



ri2i 
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Valence quark mass dependence of RRn 
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Figure 5: Numerical results of the qq+q Collaboration on 16 3 ■ 32 lattice 
at (j3 = 5.1, k = 0.177/- one parameter fit of (RRn — 1) = Xs(l ~~ £ + 
logO/(327r 2 ) ("pure chiral log"). 

with some (known) constants c k . The corresponding relations for the coefficients Wsk 
are: 2 

W Sk = W k -d k log(xs) = W' k -d k log(4xs) • (122) 

Note that these formulas can be extended to the NNLO order, too. 

A first comparison of these formulas with numerical Monte Carlo results has been 
performed by the qq+q Collaboration (DESY-Miinster) [35]. The lattice sizes were 16 
and 16 3 • 32, and N s = 2 light quark flavours were simulated. The lattice spacing 
was: a = 0.189(5) fm ~ (1.04GeV) _1 giving lattice extensions L ~ 3fm. The pion 
masses were: am w = 0.6747(14), 0.6211(22), 0.4354(68), 0.3676(23) which correspond 
in physical units to m n ~ 702, 646, 452, 415 MeV. The sea quark masses were 
approximately 60 MeV to 25 MeV; and the valence quark masses: \m sea < m va i ence < 
2m sea . Being the first exploratory study, the parameters did not correspond to the 
latest best ones, in particular, the lattice spacing was rather coarse and the quark 
masses not small enough. 

The result of this first study was that the formulas like (|118p - (|120p describe well the 
dependence on both sea and valence quark masses, in particular if some generic NNLO 
terms are included. As an example of the fits see Figure [5l First crude estimates of 
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L-G constants, renormalized at scale fo^/XR, gave with \R = 33.5(2.4) 

L m = 3.00(19) • 10~ 3 , (2L R8 - L R5 ) = -6.25(52) • 10~ 4 . (123) 

From the sea quark mass dependence it was obtained 

(2L m + Lbb) = 4.34(28) • lO" 3 , 
(4L m + 2L m - 2L m - L R5 ) = -9.1(6.4) • 10~ 5 , 

^ = 6.51(57), 
Jo 

-± = 22.9(1.5) (124) 
Jo 

These numnbers can only be taken as crude estimates, because they come from a point 
with coarse lattice spacing and no continuum extrapolation has been performed. 

5 Outlook 

The present goal of numerical Monte Carlo investigations is to perform dynamical quark 
simulations with light quarks in large volumes. After about twenty-thirty years of hard 
work - which can be considered as the preparation - the presently available computer 
resources and algorithmic developments make this goal achievable. The big question is, 
can we validate QCD as the true theory of strong interactions by comparing the results 
with experimental knowledge? After this will be done, lattice gauge theorists will be 
able to extend their research area to study at the non-perturbative level a broader class 
of Quantum Field Theories not just QCD. 

5.1 Beyond QCD 

The further development of lattice regularized Quantum Field Theories will reflect how 
the two basic theoretical problems of the Electroweak Standard Model will be resolved 
in a "beyond the Standard Model" framework. These two problems are: 

• The triviality of the Higgs- Yukawa sector: as a consequence of appearance of 
Landau-Pomeranchuk poles there are cut-off dependent upper bounds on the 
Higgs- and Yukawa-couplings, which tend to zero for infinite cut-off (i.e. zero 
lattice spacing). 

• It is very difficult to define chiral gauge theories in lattice regularization - although 
they are required for the electroweak sector. Mirror fermion states with opposite 
chirality appear and it is difficult to separate the mass scale of the mirror fermion 
sector from the known chiral sector [36J. By including the mirror fermion sector 
the theory becomes vector-like (non-chiral) . 
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These problems become acute at the TeV scale and need some solution in a near 
future - in particular based on the experimental input expected from LHC. There are 
several ways how these problems could perhaps be solved: 

1. Supersymmetric extensions of the Standard Model: the improvement of the diver- 
gence structure due to supersymmetry (the solution of the "hierarchy problem" 
because of the absence of quadratic divergences) may solve both of the above 
problems. The mirror states could perhaps be shifted to the grand unification 
scale. 

2. Technicolor-type models based on some appropriate generalization of QCD may 
produce the low-energy chiral spectrum as bound states. The mirror fermions 
could be at the technicolor scale. 

3. Beyond QFTmodels where more dimensions beyond four appear and/or quantum 
gravity effects play an important role already near the TeV scale. 

Which one (if any) of these ways is realized in Nature is a very exciting question 
and will hopefully become clear in the not very far future. If possibility 1. is realized 
then lattice field theorists will have to work more on (at least partly) supersymmetric 
non-perturbative regularization schemes. The case of possibility 2. seems to be a 
more or less straightforward generalization of QCD. In case of 3. one probably has to 
abandon the traditional QFT framework and look for radically new approaches. 
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